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ABSTRACT 

Perpendicular shocks are shown to be rapid particle accelerators that perforin optimally when the 
ratio u s of the shock speed to the particle speed roughly equals the ratio 1/rj of the scattering rate 
to the gyro frequency. We use analytical methods and Monte-Carlo simulations to solve the kinetic 
equation that governs the anisotropy generated at these shocks, and find, for iju s ~ 1, that the spectral 
index softens by unity and the acceleration time increases by a factor of two compared to the standard 
result of diffusive shock acceleration theory. These results provide a theoretical basis for the thirty- 
year-old conjecture that a supernova exploding into the wind of a Wolf-Rayet star may accelerate 
protons to an energy exceeding 10 15 eV. 

Subject headings: acceleration of particles — cosmic rays — shock waves — supernovae: general 


1. INTRODUCTION 

The “knee” at ~ 5 x 10 15 eV in the energy spectrum 
of cosmic rays arriving at Earth defines a characteristic 
energy that can be used to constrain the physics of the 
acceleration and propagation of these particles. Ideally, a 
theory of cosmic-ray acceleration would provide a natural 
explanation of this feature. Diffusive shock acceleration 
(DSA) in supernova remnants (SNR), together with the 
associated amplification of the ambient magnetic field, is 
in many respects a convincing theory, but it fails to de¬ 
liver a simple prediction of the energy of the knee. The 
maximum energy to which it can operate lies somewhat 
below 10 15 eV, and is thought to be governed by the level 
at which the nonlinear process of magnetic field ampli¬ 
fication saturates. It is not known how higher energy 
particles could be produced by within this theory (for a 
recent review, see fBeHl20Ti ' 

Th e basic problem w as iden t ified long ago 
(|Lagage fc Cesarskvl Il983al lbl: iHillasl 120051 1: DSA is 
too slow, because particles that cross and recross a 
shock front are accelerated at a rate t ~^ c , that is second 
order in the small parameter e = u s /v: t~^ c ~ e 2 w g , 
where u s is the shock velocity, v the particle speed, 
and w g = Z\e\Bc/E is the (angular) gyro-frequency 
of a cosmic ray of charge Ze and energy E in a mag¬ 
netic field B. Given this timescale, the only way for 
a SNR to accelerate particles up to the knee is by 
generating a magnetic field at the shock front that is 
substantially stronger than that in the surrounding 
medium. Magnetic field amplification in SNR is con- 
sistcn t with the observation of thin X-ra y filaments in 
SNR (|Volk et al .1120051 : IParizot et al .1120061 ) . and emerges 
from numerical simulations of para llel shocks us in g both 
particle-in-cell and hyb ri d codes dReyi]Ie_fcJ3el ]| 2012; 
Matsumoto et al. 120131 iCaprioli fe Snitkovskvl 12014 
Matsumoto et al.l2015|h However, it appears to saturate 
at a level that does not permit acceleration to energies 
above the knee dBcll et al.ll2013ft. 

In a seminal paper. I.Iokipiil (|1987h (see also lOstrowskH 
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(IMS)) noted that this problem is characteristic of par¬ 
allel shocks, at which particles cross and and recross 
the shock by diffusing along field lines, which leads to 
a relatively long cycle time compared to the gyro pe¬ 
riod, particularly if the scattering rate is low. At a per¬ 
pendicular shock, on the other hand, the cycle time is 
close to the gyro period. If DSA were to remain valid, 
this would give an acceleration rate t~^ c ~ ew g , which 
is fast enough to achieve energies well above the knee 
(jBiermannl 11993ft . The existence of a characteristic en¬ 
ergy at 5 x 10 15 eV could then plausibly be associated 
with a change in the confinement properties of the SNR 
(lDrurvll2011t iMalkov et al.l 120131) . since the gyro radius 
of a proton of this energy in the interstellar magnetic 
field is comparable with the radius of a SNR entering its 
Sedov phase of expansion. 

Early attempts to resolve this question using a Monte- 
Carlo approach for nonrelat ivistic shocks were f orced to 
make strong simplifications. iBaring et al.l (11993ft . for ex¬ 
ample, used a guiding center approximation that implic¬ 
itly assumes the distribution is almost i ndep endent of 
gyro-phase. In later work (lEllison et al.l[l995 ) this was 
lifted, and agreement with the predictions of DSA over a 
limited range of shock obliquities was found. However, a 
large-angle scattering algorithm was used that strongly 
affects the ab ilit y to r esolve anisotropies at the shock. 
iTakahara fe Terasawal (11991ft presented a method that 
takes full account of anisotropy and non-conservation of 
the first adiabatic invariant (magnetic moment), but ne¬ 
glected the transport of particles across field lines - 
a crucia l ingredient for th e treat ment of perpendicular 
shocks. iMeli fc Biermaniil (120061) investigated both the 
spectrum and acceleration timescale for highly oblique 
shocks, and found that both quantities agreed with the 
DSA predictions even for low scattering rates. However, 
they also used a guiding center approach, which can be 
justified only if the scattering rate is high. 

Later work has lifted these artificial restrictions, al¬ 
though the regime of nonrelativistic shock speeds and low 
scattering rates remains challenging. The softening of the 
accelerated particle spectrum at a perpendicular shocks 
of speed above u s = 0.1 c was investigated quantitatively 
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by iSummerlin fc Barind (I2012T) , who used Monte-Carlo 
simulations. Using a finite-difference method to solve the 
equations that result from an expansi on of the d is tribu - 
tion function in spherical harmonics, iBell et ahl 1)2018( 1 
also noted this effect for perpendicular shocks of speed 
above u s = 0.03c. In both cases, the results are in quali¬ 
tative agreement with those presented below. 

An alternat ive _ approach, pioneered by 

iDecker k, Vlahod (11986H a nd pursued, for example , 
bv lOstrowskil (|1993d and iGiacalone fc .Tokipiil (| 1996( 1 
consists of specifying the turbulent field around a 
shock front explicitly, and examining the statistical 
properties of a large number of trajectories computed 
by numerically integrating the equations of motion. 
Evidence of an enhanced acceleration rate at oblique 
and perpendicular shocks was indeed found using this 
technique. Its advantages are that it allows insight to be 
gained at the microscopic level by studying individual 
trajectories, and has the potentia l to i nclude the effects 
of anomalous transport (lie Roux et alj|2010|) . which are 
not easily accessible in a Monte-Carlo approach. 

Our main goal here is a quantitative resolution of these 
issues. We adopt a plausible microscopic model of the 
scattering process and solve for the full angular depen¬ 
dence of the distribution function, using both a Monte- 
Carlo approach and an analytic approximation obtained 
in the limit of small shock speed and low scattering rate. 
Our main results consist of the the spectral index and the 
acceleration rate as functions of the shock speed, and the 
scattering rate. 

In Section [2] we introduce the idealized model that is 
used to describe particle transport, and, in this context, 
discuss DSA at oblique and perpendicular shocks in more 
detail. The analytic approximation is presented in Sec¬ 
tion [3] and the Monte-Carlo method described in Sec¬ 
tion [4] The main results are presented in Section 0 and 
their implications for the acceleration of cosmic rays are 
discussed in Section [6] Section [7] summarizes our conclu¬ 
sions. 


2. THE TRANSPORT MODEL 

We consider the idealized situation of a plane shock 
front propagating at constant speed into a uniformly 
magnetized plasma. Cosmic rays are energetic charged 
particles, whose gyro radius is assumed to be large com¬ 
pared to the thickness of the shock front, which can, 
therefore, be treated as a discontinuity in the plasma 
velocity. The assumptions of constant, uniform fluid ve¬ 
locity and magnetic field are justified because we con¬ 
centrate on the highest energy particles, whose energy 
density is much smaller than that of the bulk of the cos¬ 
mic rays, which, in turn, is at most comparable to the 
ram pressure of the plasma flowing into the shock. 

In general, the cosmic-ray distribution function 
f(t,x,p) is a function of time t , position x and momen¬ 
tum p. It is usual to assume that particles are continu¬ 
ally deflected by magnetic fluctuations whose effect can 
be modeled as isotropic diffusion in the direction of mo¬ 
tion, i.e., as diffusion on the sphere of the end-point of 
the unit vector p/p, in addition to gyrating about the 
ambient magnetic field. This can be describ ed by the 
Fokker-Planck equation, fe.g.. lBell et al.ll2011f) which, in 
the case of a homogeneous magnetic field B embedded 


in a background plasma that is at rest, reduces to 


df _ df 

m +v ^ f+ “‘aj 


^coll 

r i 9 

fsin 9 df ) 

1 ® 2/ l 

2 

sin 0 06 

[ s de J 

sin 2 6 df 2 \ 
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where v = p/E is the CR velocity, 9 and <f> are the spheri¬ 
cal polar coordinates in momentum space with axis along 
B 1 and i/coii is a measure of the amplitude of the fluctu¬ 
ations, which is usually called the “collision frequency”. 
The assumption of isotropic diffusion in angle implies 
that r'coii is independent of 8 and </>. In the following, 
we consider a situation in which the shock front is a dis¬ 
continuity that separates two half-spaces (upstream and 
downstream) in each of which equation dTj) holds. Cos¬ 
mic rays are assumed to cross this discontinuity without 
deflection, so that Liouville’s theorem can be used to re¬ 
late the distribution immediately upstream (at t = t + , 
x = a?+) to that immediately downstream (at t = f_, 
x = X -): 

f {t + ,x + ,9 + ,<t> + ,p + ) = f (f_,£_,6U,</_,p_) (2) 


where p ±, 9 ±, and <f>± are the upstream and downstream 
momentum coordinates that label the same momentum 
vector just as t± and x± label the same space-time point. 
These quantities are connected by a Lorentz boost, in¬ 
cluding, in general, a rotation to take account of the 
change in the direction of B across the shock. 

Equations © and © suffice to describe the particle ac¬ 
celeration process: solving the kinetic equation © in the 
presence of a moving boundary (the shock front) yields 
the particle residence times and escape probabilities, and 
connecting the upstream and downstream solutions us¬ 
ing Liouville’s theorem © implements the boost in the 
energy measured in the local fluid frame which each par¬ 
ticle experiences when it crosses the shock front. It is 
interesti ng to note that a very similar system was ana¬ 
lyzed bv lSchatzmanl (]l963f ), who, however, restricted the 
scattering to changes in phase, and considered only those 
particles whose trajectories are almost tangential to the 
shock. This limited the range of validity of the treat¬ 
ment to ryiis 1 (in the notation used below), thereby 
delaying the discovery of DSA by fifteen years. 


2.1. The diffusion approximation 

If the cosmic ray distribution is almost isotropic, the 
angular dependence in equation © can be eliminated by 
expanding the momentum dependence of f in spherical 
harmonics (j.Iokiniilll971( iKingham fc Bellll2004(l . DSA is 
based on solving the resulting spatial diffusion equation 
in the presence o f a shoc k; a p articularly useful introduc¬ 
tion is given by iDrurvI (I1983T) , and this section repro¬ 
duces the relevant results in order to facilitate the subse¬ 
quent discussion. Many analytic solutions are available, 
of which two are of particular interest here. The first is 
for a steady-state particle distribution, with no cosmic 
rays far upstream and no source term above momentum 
Pq. In this case, the distribution downstream does not 
depend on position, and is a power-law in momentum: 

f(p) oc H(p-p 0 )p~ s (3) 

s = 3 r/ (r — 1) (4) 
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where r is the compression ratio of the shock and H(x ) 
is the Heaviside (or “step”) function. Importantly, this 
solution is independent of B, so that, within the diffu¬ 
sion approximation, shocks of all obliquities, including 
exactly parallel and perpendicular shocks, produce the 
same spectral index. The second solution gives the mean 
time (t) taken for a particle to be accelerated from mo¬ 
mentum po to momentum pi at an oblique shock (see 
IDnnviriflM Eq 3.31): 


on the acceleration rate: 

C <*B (10) 

u 2 _ r — 1 

=Ws ^K7TT) 

~ O (e 2 ) w g 

commonly referred to as the Bohm limit. On the other 
hand, at a perpendicular shock (\H± = 7 t/2, £?_ = rB + ) 


(t) 


u+ — U— 


f Pl dp 
Ipo P 


«± 

U+ 



(5) 


Here, (it_) is the component of the plasma veloc¬ 
ity upstream (downstream) along the shock normal, in 
a frame in which the shock is at rest and k± are the 
corresponding components of the diffusion tensor: 


t 


-l 

acc 



1 + rj 2 1 + r 
V 2 


( 11 ) 


an d the accelera tion rate rises linearly with p when p 
1 dJokipiil Il987f) . This behavior applies for all shocks 
with cos<f>_ < 1/p, which are sometimes called quasi¬ 
perpendicular, but, for simplicity, we restrict ourselves in 
the following to exactly perpendicular shocks. 


k± = K || ± c °s 2 T k_l± sin 2 T±, (6) 


with the angle between the magnetic field and the 
shock normal. In this notation, the shock velocity u s = 
it+ and the compression ratio r = u+/u-. 

Applying the diffusion approximation to equation dTJ) , 
iJokfpiil (|1987f ) showed that the diffusion coefficients par¬ 
allel and perpendicular to the magnetic field are related 
to the collision frequency by 


TgV 

k ii = 7?_ 3 _ ’ K± 


P fgV 
1 + p 2 3 


(7) 


where r g = v/ui s and the dimensionless parameter 

p = Wg/i/ C oii (8) 


describes the degree to which the particles are magne¬ 
tized. 

If p is chosen to be a constant (independent not only 
of 6 and <f>, but also of p), eq (|Sjl, results in a mean 
acceleration rate proportional to p -1 . Discussions of 
DSA conventionally adopt this scaling together with the 
additional restriction p± > 1. Although the applica¬ 
bility of these assumptions is disputed even withi n the 
framework of the di ffusion approximation (|Shalchill2009l 
iFcrrand et all [20l4lf . we nevertheless adopt them here, 
since our discussion focuses on the validity of the diffu¬ 
sion approximation itself. 

The acceleration rate t~^ c for a relativistic particle, as¬ 
suming, for simplicity, that p is the same in the upstream 
and downstream plasmas, and that p po, is given by 
equation ([5]): 


3. APPROXIMATE ANALYTIC SOLUTIONS 

Using “mixed” coordinates, in which p (and v, now 
expressed in units of c) are measured in the fluid rest 
frame, but x and t are coordinates in a frame in which 
the shock is at rest, the transport equation © becomes: 

. df df 

(l- VzU ) — + (v z - u ) c — = 

d(j) 2 

0 

Here, the shock normal is along the 2 -axis, and we ex¬ 
press both the component v z of the particle speed in this 
direction (measured in the fluid rest frame) as well as 
the speed u of the fluid (assumed to be directed along 
the shock normal) in units of c. The flow Lorentz factor 
is T = 1/Vl — u 2 , and spatial variations in the plane of 
the shock are ignored, i.e., df/dx = df/dy = 0. The 
magnetic field is in the y-z plane, and p = cos 9, so that 
v z = v sin 9 sin 0 is a function of p and (f>. For cosmic 
rays, the particle velocity, v, is close to unity. Note that 
equation 113 is valid only for perpendicular shocks; the 
corresponding equation for subluminal shocks is given by 
iKirk &; Heavens! (I1989D . 

Solutions that are stationary in the shock rest frame 
can be found by separating the variables in either the 
upstream or the downstream region: 

f(z,p) = F(p) aie KiZ ^ ,Tc Qi(p, f>) (13) 

i 

where F is an arbitrary function of p, and the eigenvalues 
A i and eigenfunctions Qi obey 


fdf_ , _L 

T 1 dcj) 2rj 


JLh - 2\ , 1 

dp [ ^ > dp i-p 2 


*acc ={ty 
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Ug c 2 
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rjr 

cos 2 + 


2 sin T + 

cos 2 '&+ + 


1 + rj 2 


sin 2 T- 


1 + p 2 


(9) 


At a parallel shock (^± = 0, = B + ) equation ©, 

together with the restriction p > 1 gives an upper limit 


A i (v~ - u)Qi = 

d . 2] d Id 2 ' 

dp ^ A< ' dp + i - p 2 df! 2 

(14) 

together with boundary conditions on Qi that ensure reg¬ 
ularity and single-valuedness on the sphere. In general, 
equation has an infinite number of both positive 
and negative discrete eigenvalues (labeled with i > 0 and 
i < 0 respectively), in addition to the eigenvalue A 0 = 0 


j_d_ J_ 
1 d(f> 2p 
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with eigenfunction Q o = constant. These govern the 
spatial dependence of the solution: the isotropic part is 
independent of z, whereas the eigenfunctions with i < 0 
decay exponentially towards positive 2 (i.e., upstream), 
and grow exponentially downstream on a length scale 
that decreases as |z| increas es. Fo llowin g the proce dure 
used for relativistic shocks dKirk fc Schneideilll987fl the 
solutions in the upstream and downstream regions can be 
matched at the shock to find the function F(p), which, 
in the absence of a source term, is a scale-free power 
law F oc p~ s ■ The boundary condition far upstream (at 
z —> 00 for u+ > 0 ) is enforced by restricting the ex¬ 
pansion to i < 0; that far downstream is imposed by 
requiring the projection onto downstream eigenfunctions 
with i > 0 to vanish, thus preventing divergence of the 
distribution as z — > —00. 

In the case of parallel, relativistic shocks, the expan¬ 
sion m was found to converge rapidly, and an accurate 
result was obtained by retaining only a single eigenfunc¬ 
tion dKirk et al.l 12000!) . The current problem is more 
complicated, since the eigenfunctions are functions of 
gyro-phase <j> as well as pitch angle 9 , and the eigenvalue 
problem changes from a single-parameter (u s ) problem 
into a two-parameter (u B , if) problem. This makes an 
expansion to high order cumbersome. Also, the differ¬ 
ential operator in Q is not self-adjoint, so that the 
adjoint operator (obtained by changing the sign of the 
d/d<t> term) must be used to find the function needed for 
projection onto the “forbidden” downstream eigenfunc¬ 
tion. Nevertheless, the analogy is close, so that one can 
again hope to find a reasonable approximation by using 
only a single eigenfunction. Adopting this approach, the 
power-law index s for a shock moving at speed it+ into 
the upstream fluid and at speed u_ in the downstream 
fluid is given implicitly by the equation 

S = JJd^.Q. (/z_,0_) (v z _ — U-) x 

(p+/p-)~ s Q+ (P+A+) 

= 0 (15) 

where the label i = — 1 of the retained, leading eigen¬ 
function has been omitted, and the notation Q is used 
to denote the corresponding eigenfunction of the adjoint 
equation. The notation (n±, <f>±) is used to indicate an¬ 
gles in the up and downstream fluid frames. For particles 
that move rapidly with respect to the shock, m can be 
expanded to first order in u/v, leading to a linear equa¬ 
tion for s, with solution: 


two respects. Firstly, since it falls off more slowly with 
distance from the shock than do the other eigenfunc¬ 
tions, it is the dominant contribution to the distribu¬ 
tion function far upstream, and so can have no zeroes 
in the range 0 < (f> < 27 t , — 1 < fi < 1 . Secondly, in 
the limit u — > 0 , it merges with the eigenfunction i = 0 , 
which is a well-known feature characteristic of the dif¬ 
fusion approximation dFisch fc Kruskal il980f) . Taking 
this limit, it is straightforward to show that both A_i 
and the anisotropic part of Q-\ are of first order in u. 
The transformation of the arguments of the eigenfunc¬ 
tion Q + in equation (flfil) from n + ,<j> + to acts on 

the anisotropic part of this function only. In the case 
of highly relativistic particles, to which we restrict our¬ 
selves in the following, it produces a modification of this 
term that is of first-order in (u+ — u_). It follows that 

Q+{p+A+)=Q+(p-,<t>-) + 0 (u+ 2 ) (17) 

and, using the orthogonality relations 

jj^ + (v z+ -u + )Q_=0 (18) 

dp-dtp- (v z _ — U-) Q+ = 0 (19) 

in equation (fl 6 l) , leads immediately to the standard DSA 
result ( 0 . 

However, this analysis is based on the assumption that 
u is the only small parameter in the problem, which 
breaks down when scattering is sufficiently weak. As¬ 
suming the ordering u ~ 1 / 77 ~ e <C 1 , standard pertur¬ 
bation techniques (for details see appendix [A} lead to an 
expression for Q that is anisotropic at zeroth order: 

Q = a(p)e Av ^ T -^ cos ^ + O(e) ( 20 ) 

where 

a(p) = Ps°( M ,-A 2 /2) ( 21 ) 

with Psj] the angular, oblate, spheroi dal wave fu nction 
with m = n = 0 (in the notation of iThompsonI I 20 TTI . 
chapter 30). The eigenvalue A is related to the spheroidal 
eigenvalue A(J by 

Ag (-A 2 /2) = A (A + 2 77 it) (22) 

and is of zeroth order in e. 

Evaluation of the <fi integrals in (fl 6 l) is straightforward. 
The integrals over p, do not appear to be possible ana¬ 
lytically, but are simple numerical quadratures. An ex¬ 
pansion in powers of r/u yields 



' ” («+ - «-) X 

// dp_d(/)~Q~ (/r-,<M (vz- ~ u~)Q+ (P+A+) 
If dp,-d<p-Q- ( p-,(/)-)vj_Q + (a*+,<£+) 

(16) 

where, for simplicity, the relative speed of the down¬ 
stream fluid with respect to the upstream fluid is 
assumed to be u+ — rt_, implying that they both 
flow along the shock normal, which is the case for a 
shock of high Alfveni c mach number (see, for example, 
iKirk fc Heavensl[l989f) . 

The retained eigenfunction with z = — 1 is special in 



9 (r + 1) 
20r(r — 1) 


7 f u + 2 + O (7? 4 u + 4 ) 


(23) 


indicating a significant softening of the spectrum for fi¬ 
nite r]u+, as compared to the standard DSA result. In 
section [5] we compare this result and the angular depen¬ 
dence of the eigenfunction defined by equations (l 20 l) and 
(l2ll) with the results of the Monte-Carlo simulations de¬ 
scribed below. 


4. MONTE-CARLO SIMULATIONS 

To solve the fully relativistic kinetic equation ©> 
we use two conceptually different, but closely related 
Monte-Carlo methods, similar to those described by 
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lAchterberg et al.1 (1200111 and bvlEllison fc Doubld (|2004f ) 
and iSummerlin &; Barind (|2012iK 

In the first, we rewrite m in the form of a Fokker - 
Planck equation, and use a theorem due to lltoh ( Il944f l 
to write down a stochastic differential equation governing 
a family of effective trajectories, whose statist ical prop¬ 
erties are those of the required solution f (see iGardineil 
fl99l chapter 4). An explicit first-order discretization 
scheme is then used to construct a large number of these 
trajectories. In the second, we start from the Boltzmann 
equation describing a distribution of particles that move 
in the unperturbed, uniform magnetic field, but are sub¬ 
ject to random collisions, each of which causes a small 
but finite angular deflection. In the limit of small and 
frequent deflections this equation can be reduced to ©, 
as we demonstrate in Appendix [C] but we solve it for 
small finite deflections, advancing the trajectory in be¬ 
tween the scatterings by numerical integration. Details 
are provided in Appendices |B] and [C] We have verified 
that the results presented in Section [5] do not depend on 
which algorithm is employed. 

In each case, trajectories are initiated at the shock 
front. They then perform an excursion that either re¬ 
turns it to the shock, or terminates the trajectory when 
it crosses a boundary placed at a large, fixed (in units 
of the gyro-radius r g ) distance from the shock front in 
the downstream region. An excursion that returns to the 
shock front initiates a subsequent excursion in the other 
half-space, starting at the same space-time point with the 
same space-time coordinates and momentum four vector. 

The anisotropy of the particle distribution at the shock 
front is the characteristic feature of this problem. In ad¬ 
dition, we are interested in two properties of the solu¬ 
tions: the time-asymptotic power-law index at momenta 
well above injection, and the mean time taken for ac¬ 
celeration to a given momentum. These quantities can 
be used to determine whether or not a SNR lives long 
enough to enable acceleration at quasi-perpendicular 
shocks to play an important role. Each of them can be 
extracted from a simulation of a large number of trajec¬ 
tories, all injected at momentum po (assumed me) at 
time t = 0, by recording the sets of values ti and p) at 
each crossing of the shock front. These values, collec¬ 
tively called events , are labeled by the integer i. 

The distribution at the shock front as a function of the 
angles 9 and <j), and the momentum p (in each case inte¬ 
grated over the other two variables) is obtained by set¬ 
ting up logarithmically spaced bins, and adding to these 
a weighting factor equal to the reciprocal of the relative 
velocity of the particle with respect to the shock surface. 
Particles are initiated at the shock with an angular dis¬ 
tribution that is arbitrarily chosen to be isotropic in the 
hemisphere entering the upstream plasma. However, af¬ 
ter a few shock crossings, the bias thereby introduced is 
lost, and the distribution settles down to one that is inde¬ 
pendent of the number of crossings, and can be compared 
to the analytical form found for the scale-free, power-law 
distribution. Provided the average number of crossings 
per trajectory is large (in the example presented below 
it is roughly 100) the effect of the bias is not noticeable. 
After a few crossings, particles display a power spectrum 
which extends from slightly above po U P to a momentum 
determined by the position of the downstream boundary, 
or the time limit placed on the trajectory. Within this 


range, the time-asymptotic power-law index s is found 
from a least-squares fit in log-log space. To find the av¬ 
erage acceleration time, the events are again binned in 
p, and, at the same time, a running average of ti is ac¬ 
cumulated. In the results presented below, s and ti were 
determined using the range 0.1 < log 10 (p/po) < 1- 

5. RESULTS 

We present the results of fully relativistic simulations 
of perpendicular shocks for an upstream speed u s (= rt + ) 
ranging from 0.01 to 0.2 (in units of c), and a turbulence 
parameter 77 between 1 and 100, and compare them to 
analytic approximations. A total of 50,000 trajectories 
were simulated for each set of parameters. For conve¬ 
nience, the plasma is assumed to flow along the shock 
normal both upstream and downstream, and the com¬ 
pression ratio r = u+/u- is chosen to equal 4, which 
corresponds approximately to the value derived from the 
Rankine-Hugoniot conditions for an ideal gas of specific 
heat ratio 5/3, (although it is not relativistically exact for 
any physically motivated equation of state). The mag¬ 
netic field strengths B+ and R_, measured in the up 
and downstream rest frames, respectively, are related by 
B + /B- = rT s ^/l — ( u s /r ) 2 (see iKirk fc Heaveiislll989L 
eq (4)). Itoh’s method was used in the simulations pre¬ 
sented here. 



(j)/-K 



Figure 1. The angular dependence of the distribution function / 
at the shock (arbitrary normalization) for u s = 0.012 and 77 = 22 . 
Monte-Carlo simulations (red points) are compared to zeroth-order 
analytic approximations. The top panel shows the phase depen¬ 
dence, (/ integrated over pitch angle), the bottom panel the pitch 
angle dependence, (/ integrated over phase). The bl ue l ines are 
found by integrating the expression given in equation pot numer¬ 
ically. 

Figure [I] shows the angular distribution at the shock 
front found from Monte-Carlo simulation as a function 
of gyro-phase (f> (top panel) and cosine p of the pitch an¬ 
gle (bottom panel), for upstream speed u s = 0.012 and 
turbulence level rj = 22 (r]u s = 0.26). Particles are reg¬ 
istered by the simulation as they cross the shock front, 


























6 


Takamoto & Kirk 


so that very few events accumulate in bins where the 
velocity vector lies close to the plane of the shock, i.e., 
at \J\ — fi 2 sin <f> r; u s . This accounts for the relatively 
large fluctuations close to (f> = 0, ± 7 r, which corresponds 
to grazing incidence for most values of /i. The maxi¬ 
mum of the distribution function lies close to (j> = ±7r, 
which corresponds to particles moving along the shock 
front in the direction of the shock- or “grad-B”- drift. 
The zeroth-order analytic approximation reproduces the 
main features of the simulated distributions, although it 
is more strongly peaked both in <f>. This may be either be¬ 
cause the eigenfunction expansion requires several terms 
in order to converge to an accurate solution, and/or be¬ 
cause the the analytic approximation is, strictly speak¬ 
ing, valid only in the limit of nonrelativistic shocks and 
low scattering rates ( u s ~ 1/rj —> 0). In any case, it is 
clearly important to employ a scattering algorithm in the 
simulations that resolves angular structure on the small 
scale indicated by the eigenfunction. 



-1.5 -1 -0.5 0 0.5 


loglO (r/Us) 

Figure 2. The energy spectral index s (top panel) and the ac¬ 
celeration rate divided by the DSA prediction (bottom panel), as 
a functions of r]u B . Red crosses show the results that fall in the 
plotted range of Monte-Carlo simulations with 10 values of ?/. be¬ 
tween 1 and 100, and 30 values of u s between 0.01 and 0.2, in 
each case equally spaced in logarithms. Cyan and green curves in 
the top panel show ana lyti c approximations using the series expan¬ 
sion given in equation (1231) , and numerical integration respectively. 
Blue curves show fits to the simulation results. 

Figure [2] shows the time-asymptotic power-law index 
s (top panel) and the average acceleration rate (bottom 
panel) at a perpendicular shock front, as functions of 
the product tju s , for various values of the parameters u s 
and rj. In each case, the results fall close to a single 
curve with very small dispersion about it. This suggests 
that the most important parameter in this range is the 
combination r]u s , which is also the only parameter of the 
approximate solutions found in section [3] Physically, r]u s 
is the collision time divided by the time taken for the up¬ 
stream flow to travel across one particle gyro-radius, i.e., 
roughly the inverse of the number of collisions expected 


as an undisturbed particle orbit advects across the shock. 
Provided the fluid speed is nonrelativistic, this is the only 
physically significant dimensionless quantity in the prob¬ 
lem as we formulate it. 

For r]u s -C 1, the Monte-Carlo simulations show a spec¬ 
tral index s that lies close to the DSA prediction s = 4, 
in agreement with the analytic approximations. The ac¬ 
celeration rate in this regime is also close to the DSA 
prediction, which is a factor 2.5 r/ (l + l/f? 2 ) faster than 
the Bohm rate, as defined in equation m- Inspection of 
the angular distributions shows, as expected, that they 
are almost isotropic in this region. 

However, as rju s increases, simulations show that the 
index s increases (i.e., softens). This is also seen in the 
analytic approximations, although the softening here sets 
in at somewhat higher rju s , and proceeds more rapidly. 
According to the simulation results, a softening of unity, 
i.e., s = 5 is reached at r/u s ~ 1, whereas the zeroth-order 
analytic approximation predicts this index at r]u s ~ 2 . 
The softening of the spectrum is accompanied by a re¬ 
duction in the acceleration rate, expressed in units of the 
DSA-predicted rate. At r]u s = 1, the reduction is roughly 
a factor 2 , implying that acceleration still proceeds at a 
rate that is about a factor of ry faster than the Bohm 
rate. 

6 . DISCUSSION 

The results presented above demonstrate that per¬ 
pendicular shocks impose a strong anisotropy on parti¬ 
cles that are scattered in the upstream and downstream 
plasma, when the ratio 1 /77 of the scattering rate to the 
gyro frequency is less than or of the order of the ratio 
of the shock speed to the particle speed. Instead of dif¬ 
fusing in space around the shock front, as predicted by 
DSA, the accelerated particles are then concentrated into 
a fan-beam structure within a relatively small interval 
A <p of gyro-phase directed in the plane of the shock. This 
can be seen from the phase-dependence of the eigenfunc¬ 
tion given in equation (l20l) : Q oc exp[Aw-\/l — /i 2 cosc/>]. 
The direction of the beam corresponds to the drift im¬ 
posed on an unscattered trajectory by the presence of the 
shock front. For r)u s > 1, the opening angle of the fan 
beam can be estimated from the asymptotic expression 
for the eigenvalue, A ~ —3ryu s (see Appendix 0, to be 

A (j) ~ (rjUs)^ 1 ^ 2 . This roughly corresponds to the diffu¬ 
sive spread in phase produced by the scattering operator 
in equation ©, when acting on a perfectly collimated 
beam over the time (« 1/ ( uj g u s )) taken for an unper¬ 
turbed trajectory to cross the shock. 

In this regime, a particle that is moving close to the 
shock front on its downstream side which finds itself in¬ 
side the fan beam has a relatively high probability of 
recrossing the shock front many times. On the other 
hand, a particle in the same position but with a momen¬ 
tum directed out of the beam is likely to be swept away 
downstream after only a few crossings. For this reason, 
the intuitive picture of DSA, which is based on assign¬ 
ing all accelerated particles in the downstream plasma an 
escape probability that is independent of their position 
and direction of motion, is inadequate. The more for¬ 
mal derivation of DSA based on the diffusion-advection 
equation also breaks down, because spatial diffusion de¬ 
scribes the transport process only when the distribution 
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is approximately isotropic. 

We find that the gradual breakdown of DSA at a non- 
relativisic, perpendicular shock as the collision rate de¬ 
creases depends on the single parameter rju s and has two 
important effects on the accelerated particles. Firstly, 
it softens the spectrum. For a compression ratio of 4, 
the phase-space density, / oc p~ s , has s ~ 0.7r]u s + 4, 
compared to the DSA prediction of s = 4. Secondly, 
it causes a reduction in the acceleration rate by a fac¬ 
tor of approximately l.l? 7 u s + 1 compared to the DSA 
prediction, which rises linearly with 77 . Thus, optimal 
conditions for acceleration, in the sense that it proceeds 
rapidly and produces a reasonably hard particle spec¬ 
trum, are found when the anisotropies induced by the 
shock speed and the magnetic compression are compara¬ 
ble, i.e., when iju s ~ 1. 

The validity of DSA at perpendicu lar sh ocks has been 
a controversial issue for many years. I.Tokiniil (Il987h sug¬ 
gested that breakdown would happen when the collision 
time becomes longer than the time during which a gy¬ 
rating particle interacts with the shock, leading to the 
approximate condition rj < l/u s and an upper limit on 
the acceleration rate that is first order in e: 

Me < V Mb) 

r — 1 

= U s - CUa- 

r 

~ O (e) u g . (24) 

This estimate is in rough agreement with our findings. 
lAchterberg fc Ball (1994), on the other hand, proposed 
that the collision frequency must be large enough to al¬ 
low particles to diffuse along the magnetic field whilst 
upstream of the shock. This leads to the restriction 
rj < 1 / y/ul and 

M < M /2 M 

= u s /2L y^ U} g 

~ O (M) (25) 

Our result that the acceleration time and spectral index 
are functions of the product r)u s does not support this 
conjecture. However, the situation may be different at 
oblique shocks, where particles have more opportunity to 
diffuse along the upstream magnetic field lines. 

The importance of these effects for the acceleration of 
high en ergy cosmic ra ys has b een emphasized i n partic¬ 
ular by iJokipiil (119871) and by iBiermannl (11993H . Never¬ 
theless, cosmic-ray acceleration at quasi-perpendicular 
shocks has received relatively little attention for several 
reasons. Firstly, observations of the polarization of radio 
emission in SNR 1006 show that accelerated electrons 
are found predominantly in regions where the magnetic 
field is disordered, and where the normal to the shock 
front lies roughly in the direction of the external field , 
rather than perpendicular to it (jRevnoso et al.l [2013D . 
This general trend is also consistent with the predomi¬ 
nantly radial orientation of the i nterior magnetic field in 
young supernova remnants dRev nolds fe Gilmorel 119931 : 
iParizot et al.1120061 iRevnolds et al.l I2012D . and fits in 
with the idea that magnetic field generation acts at 
quasi-parallel shocks. Secondly, numerical simulations 


of ac ce leration using Mo n te-Ca rlo codes (lEUisoiL_et_aI 


1995; lEllison fc Double! 120041 iSumm erlin fc Baring 
2012 ). hybrid co des (TCaprioli fc Spitkovskvl 2014T) 


and PIC codes (IStockem et al.l I2012~E ISironi et al.1 
120131 IMatsumoto et al.l 120 1 31) also disfavor the quasi¬ 
perpendicular orientation, which, according to these 
results, is less efficient at injecting particles into the ac¬ 
celeration process, provided the initial level of turbulence 
is very low. 

However, these reasons do not directly apply to the 
problem of the acceleration of very high energy ions, on 
which we focus our attention in this paper. On the one 
hand, observations of synchrotron radiation relate exclu¬ 
sively to the electron distribution, which is likely to be 
much more tightly confined to the shock front than are 
high-energy ions. On the other, numerical simulations 
consider an initially uniform upstream magnetic field. 
This may indeed inhibit injection, but, in a more realis¬ 
tic situation, some level of turbulence must be present 
initially. Our results indicate that efficient accelera¬ 
tion can be expected for an effective collision frequency 
z'coii ~ u s ujg, which implies a very low level of turbu¬ 
lence under SNR conditions. However, if the collision fre¬ 
quency is even lower, we find that perpendicular shocks 
should produce only very steep spectra, in agreement 
with the apparent failure of these shocks to trigger ac¬ 
celeration in numerical simulations. Furthermore, in a 
realistic situation, low energy particles may see localized 
regions of quasi-parallel geometry due to small length- 
scale fluctuations in the upstream medium and/or the 
shock speed. These regions may be very effective accel¬ 
erators, possibly giving rise to substant ial field amplifi¬ 
cation. But, as described bv iBelll (12014H the affected ac¬ 
celeration region remains limited in spatial extent, and 
automatically provides an escaping flux of particles at 
the highest energy to which it operates. In the context 
of the problem we consider here, these can be regarded as 
being injected into an acceleration process operating on 
larger spatial scales, on which the shock is perpendicular. 

Another well-known argument against acceleration to 
energies above the knee at perpendicular shocks is that 
particles drift across the shock surface whilst undergoing 
acceleration, covering a distance proportional to their en¬ 
ergy. Depending on the geometry of the field, this might 
move them out of the region where the shock is perpen¬ 
dicular. If one sets an upper limit to this distance equal 
to the radius of the SNR, the corresponding upper limit 
on the energy turns out to be comparab le to that found 
for quasi-parallel shocks (e.g.. lBclll2014fl . and is roughly 
equal to the energy a particle could gain by drifting from 
the pole to the equator (or vice versa) in a flow moving at 
the shock speed through a steady, magnetized, axisym- 
metric stellar wind. 

However, whether or not particle drift really limits the 
maximum energy depends on the specific configuration 
of the magnetic field. For example, in a uniform exter¬ 
nal field the drift motion induced by a spherical shock 
front is directed along lines of constant latitude (mea¬ 
sured on the shock surface with the polar axis along the 
direction of the external magnetic field), along which the 
shock does not change its obliquity. Diffusion along the 
magnetic field lines in this configuration might allow par¬ 
ticles to escape to regions where the shock is parallel, but 
this effect depends on the transport properties in the up- 
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stream and downstream plasmas, rather than on the drift 
motion. 

The situation is different for a spherical shock front 
expanding into a magnetic field that is anchored in the 
wind of a rotating progenitor star. In the ideal case of 
a magnetic dipole aligned with the rotation axis, a small 
region with quasi-parallel configuration may exist close 
to the axis. However, if the progenitor is a miss-aligned 
rotator, and/or one that undergoes repeated field rever¬ 
sals, the undisturbed field is perpendicular to the shock 
normal essentially everywhere, once this has expanded 
to well beyond the Alfven radius. In this case, the direc¬ 
tion in which a particle drifts as the shock moves over a 
magnetic field that is frozen into the progenitor’s wind, 
can be a rapidly changing function of the shock radius. 

The majority of supernovae are thought to result from 
the explosion of massive stars that can be assumed to 
have had a strong wind, and may have been both mag¬ 
netized and rapidly rotating. Assuming a wind speed of 
u w , and an Alfven surface not too far from the star, the 
toroidal magnetic field at large radius R can be estimated 
as B(R) = B* (f2R*/i> w ) (R*/B), where I?* is the stellar 
radius, B* the surface magnetic fi eld and Q the angular 
velocity of the star (iParkef 1958 ). Then, defining the 
maximum energy B max to which a proton can be accel¬ 
erated by identifying the acceleration rate given in equa¬ 
tion ED with u s c/R , setting r = 4, an d insert i ng values 
thought typical of Wolf-Rayet stars dBerezhko fc Volki 


1200ft . leads to 


E 


max 


3 

8^ s 




»w 


eB*B* 


1.710 lf Vs 




R* \ 
3.10 12 cm / 


eV, 

(26) 


We find that equation ED over-estimates the accelera¬ 
tion rate by a factor of 2 when rju s ss 1, which is unim¬ 
portant in view of the uncertain values of Q and B*. 
Therefore, we conclude that, under optimal conditions, 
i.e., when rju s « 1, such objects may be capable of ac¬ 
celerating protons into a power-law spectrum somewhat 
softer than that predicted by DSA up to energies well in 
excess of 1 PeV. 


7. CONCLUSIONS 

We have studied the distribution function of particles 
accelerated at perpendicular shocks as their scattering 
rate decreases, finding that a strong anisotropy devel¬ 
ops, which causes the theory of diffusive shock accelera¬ 
tion to fail. When the scattering time is comparable to 
the time taken for a particle orbit to traverse the shock 
front, the power-law index of the particle distribution is 
found to soften by roughly unity, compared to the DSA 
prediction and the acceleration rate is roughly halved. 
These results provide a firm basis in kinetic theory for 
the long-standing conjecture that protons can be accel¬ 
erated to energies well above 1 PeV at the perpendicular 
shocks that are expected to form when a supernova ex¬ 
plodes into the wind of a massive progenitor. 
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APPENDIX 

APPROXIMATE ANALYTIC SOLUTION 

In equation ( 1141 ) we assume u^l/? 7 ~e-Cl, reintroduce the label i, and pose an expansion: 


Qi ->Qi ^ + tQi + 0(e 2 ) 

(Al) 

Aj — >• A- 0) + eA- 1 ) + 0(e 2 ). 

(A2) 

Since v z = VyJ 1 — /z 2 sin^>, the zeroth order terms 


dQ^ A (0) q(0) _ q 

d</> + 1 

(A3) 

can be integrated to give (writing /z = cos 6 where it makes the notation more compact): 


Qf ] = ai (/z) e A i 0> sin0cos ^ 

(A4) 

and, for the adjoint function 


Q[ 0) =di (/z) e- A ‘ 0> sin0cos ^. 

(A5) 
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The first order terms are: 


dcj) 


+ hJf' 1 sin 9 sin cf> 


- A s (1) sin 9 sin <j?J (j( 0) + 


1 


d 2 


d , 2 v d 1 

2r] [<9/x ^ djj, + 1 — /x 2 dcj) 2 


Q 


( 0 ) 


Integrating, using the integrating factor e A ( 1 Sln e cos ^ gives 


Q\ = e 


sin tf cos c 


F(cj>) = e- Ai sin6cos 

1 


r d<f/F {cj>') + (/x) 

Jo 4 ry 

f^u — A - 1 ' sin 9 sin cj?J Q + 


d 2 


d , 2 , d 

2 V \djl ' A< ' dfl + 1 - fi 2 dcj) 2 

Imposing periodic boundary conditions in <j ), on Q ^ leads to 


Ql 


(o) 


d / 2\ dcii A, 

which is to be solved with boundary conditions 

A (°) 


(o) 


ai (1 ~ p2) ij + if ( A '“ > (1 + ^ + M °*= °’ 


da* _ 

d/i ”^2 


^A- 0 ^ + 2rjuj a,i at /x 


= ± 1 . 


The first-order terms are then: 


Qf } = - 


A^ sin 6 cos <f> 


4 rj 

Aj 0) sin 9 sin cj> 


b + sin 9 (cos cj> — 1) a* — 

^4 + A ? - 0 ^ sin 9 cos tpj aj + 4 cos 0^- 




g—A^ sin 0 cos <f> 

4rj 

Aj- 0 -* sin 9 sin <j> 


b — 477 A sin 9 (cos cj) — 1) a* — 


4 — A.- 0 - 1 sin 9 cos 


where b ~ 1 and b ~ 1 are functions of /x which, together with A^ 1 , can be constrained by examining the second order 
equations. However, because these terms are even in q i>, they do not en ter into the comp utation of the index s. 

Equation (IA9I) is a special case of the spheroidal differential equation (lThompsonll201lL chapter 30), whose solutions, 
for real A, : , are the oblate, angular, spheroidal wave functions Ps™ (/x, — A?/2). In general, solutions exist for both 
positive and negative A$, together with the special isotropic solution Ao = 0, ao = constant. However, to represent 
the upstream distribution we require an eigenfunction which has no roots for — 1 < /x < 1. This is the function 
Psq (/x, — A^/2). The eigenvalues of the spheroidal wave equation A™ (q 2 ) (in the notation of lThomnsonl (I2011IH are 
defined for n > m and ordered such that A™ < A™ +1 <_They are related to A, by 


, „&CLi 

CLi + 4cosfl —— 
d ii 


(A6) 


(A7) 


(AS) 


(A9) 


(A10) 


(AH) 


(A12) 


A° i_i ( — A 2 /2) = Aj (Aj + 2rju) . 


(A13) 


Thus, as expected, the largest negative eigenvalue A_i corresponds to Aq and, therefore, to the eigenfunction that has 
no roots in the range —1 < /x < 1. For large 77 /x., A_! —> —2r]u. A power series in r/u, can be found directly from (IA9I) 
(see also lMeixner fc Schafkdfl954 l p 240): 


A = —3r]u + — (r]u) + O (rfu 5 ) 

a = 1 + ^ (vu) 2 + ( 4 + 9/^) (wf + O (rfvP) 


(A14) 

(A15) 


Using these series to evaluate the integrals in equation (flbl) leads to equation 
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MONTE-CARLO METHOD BASED ON STOCHASTIC DIFFERENTIAL EQUATIONS 
Equation 0, when rewritten in Fokker-Planck form, becomes: 

\ d d 2 d 2 ( f 

) dcf) dp 2 y 2r] J dcj) 2 \2r)(l — fj, 2 ) 

where p = cos 9, time and space are measured in units of uj ~ 1 and c/uj s , and v is in units of c.. Solutions to this 
equation can be found by constructing sample trajectories that satisfy the set of stochastic differential equations (see 
lGaidinejH99l 



df _ d d 

dt dx ^ V ft + dp 


V 


dx = vdt 

dp = - ( p/r )/) dt + [(1 - v 2 ) /r )\ 1/2 dW t (B2) 

d^ = dt + [(i - p 2 ) rj\ _1/2 dw t 

where dWt is an infinitesimal Wiener process. We use an explicit first-order discretization to find a numerical solution 
to this set, consisting, for each sample trajectory, of a sequence of points in phase space labeled by i: 


Xi + 1 = Xi + UjAt 

Pi+i = Hi - {Hi/rj) At + £,i\ 


'At (1 — /if) 


0i+i — 0i + At + 


At 


(X~ P 2 )v 


(B3) 


where and Q are random numbers uniformly distributed on the interval (—\/3, a/ 3) , which, therefore, have zero 
mean and unit variance. This scheme is rapid, but has the disadvantage that trajectories that pass very close to the 
points p = ±1 are subject to errors. The affected range depends on the time-step, and is given approximately by 
|/x| >1 — At/r). Typically, we choose At = 10~ 2 and 1 < r) < 20, so that < 1% of particles in an isotropic distribution 
are affected. This is unimportant in the simulations presented above, where no fine-scale structure in /i is expected. 
However, it could become a concern for parallel, relativistic shocks with T > 100. 

Each sample trajectory starts on the shock front and consists of a series of excursions into the up and downstream 
plasmas, ending when it crosses a boundary placed a fixed distance d downstream of the shock. We performed 
simulations for several different values of d. For oblique shocks, d ~ 200 is sufficient, in the sense that a power-law 
distribution in p is reproduced over several decades. However, parallel shocks require d ~ 1000, reflecting the fact that 
the diffusion length along the magnetic field is greater than that across it, if i) > 1. When a timestep causes a trajectory 
to cross the shock front, the step is repeated with smaller At chosen to place the particle precisely on the shock front. 
The total number of timesteps in the excursion is then checked, and the excursion repeated if this number is too small 
— typically less than 3. This procedure eliminates trajectories that depend strongly on the finite size of the timestep, 
at the expense of distorting the angular distribution of particles at the shock front that move almost tangential to it, 
thereby limiting the ability of the simulation to resolve fine-structure in gyro-phase 0 in directions close to the plane 
of the shock. From the analytic approximation, such structure should be present on the scale A0 ~ 1/A ~ (t]u s ) , 
so that the choice At = 1/100 limits the accessible parameter range to rju s < 100. 


MONTE-CARLO METHOD BASED ON THE BOLTZMANN EQUATION 


In this method we solve the equation that results when the right-hand side of equation 0 is replaced by the 
Boltzmann collision operator corresponding to elastic pitch-angle scattering in the test particle approximation, given 

by 


d£ 

at 


coll 


i 


du'dtf f(n', 4>')p{n, />; //, 0') - /(//, 0) 


-1 ^0 


(Cl) 


where £ sca tt is the mean time between scatterings and p(/x, 0; p!, 00 is the probability that a parti cle with (u/00 is 
scatte r ed into (m 0) in a singl e scattering. This equatio n is solved using the method developed bv lEllison fc Eichlerl 
(I1984T) : Ijones fe Ellisonl (I1991H : lEllison fe Doublel (I2004T) . Trajectories whose statistical properties are given by / are 
found by numerically integrating the standard relativistic equations of motion in a steady background electromagnetic 
field — corresponding to the left-hand side of equation 0 — using the exact solution or the Bulirsh-Stoer method 
over a fixed time interval chosen to equal the average time between scatterings: £ sca tt = 2-7T /w g N, where N (3> 1) is a 
parameter of the scattering model. (Formally, a random, exponentially distributed time interval with mean value f scat t 
should be employed, but this is not necessary in the present problem, where the number of scatterings between each 
encounter with the shock is required to be large). On scattering, trajectories are su bject to a random deflection t hrough 
a small angle that is uniformly distributed between zero and O max (for details see iSummcrlin fc Baring] (120121 11. This 
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enables one to identify a timescale ti so on which the particle distribution would, in the absence of a driving mechanism, 
relax to isotropy: 

tiso = 6t SC att/0 max (02) 

The equation of motion is solved in the upstream and downstream half-spaces in the respective comoving frame, 
allowing the sample trajectories to cross the shock front without deflection or energy change. The escape boundary is 
located in the downstream region where the probability of particles returning to the shock front is sufficiently small, 
as described in Appendix iBl 

In order to demonstrate the link between this method and the stochastic differential equation approach described in 
appendix [Bj we note that when the limit N —> oo is taken with ti so fixed, the scattering model corresponds to a phase 
function 

p (/q 0;//, P) =p{ cos©) 

= H (cos 0 max - cos 0) / (27r0 max 0), (C3) 

where we assume 0 <C 1 and use the normalization ff d(j)dQ sin 0p(cos 0) = 1. Expanding this phase function by 
writing 

oo n 

P( cos0)=47r£ £ a n Y™(0,p)Y™(6',p), (C4) 

n—0 m=—n 


where Y™ is a spherical harmonic, and the asterisk means the complex conjugation, one finds 


1 

<Jn ~ 2 


0 2 


12 


i{n + 1) 


Using equations (IC4I) and (IC5I) . the collision operator, equation (IC1I) . can be rewritten as follows: 

a2 r 1 /*2 tt 

^max 


d£ 

dt 


12 t s , 


dii'dp 


coll i ^<'sca,tt J—lJo 

where we used the completeness relation: 


oo n 


£ £ n(n+l)Y™(e,<f>)Yp™(e',p) 


n —0 m——n 


fit*, P) 


£ £ Y™(0, P)Y~ m (d', P) = < 5 (m - P)S(<I> - P). 

n =0 —n 

Furthermore, using the following two relations: 


— n(n + l)Y™ = 


d M 2 .d id 2 


f»l r 27T 


dndtf(p)Y™ = - 


1 


-1 Jo 


n{n + 1) J_ 1 


dfi ' 1 — p 2 dp 1 

r l c 2tt 

d/idpY™ 


Y n 


d f1 2 .d id 2 


the collision operator finally becomes: 

'dj_ 

dt 


coll ^iso 


dfx dfi 1 — p 2 dp 2 _ 

d 1 d 2 


d^ 




p 2 dp 2 


WP), 


which is equivalent to ([lj one if we identify ti so = l/^ C oii, or, equivalently, t? = u] g ti so . 


(05) 

(C6) 

(C?) 

(C8) 

(C9) 

(CIO) 
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